{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://csdms.colorado.edu/wiki/ESPIn2020\"><img style=\"float: center; width: 75%\" src=\"../../media/ESPIn.png\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diffusion\n",
    "### Introduction\n",
    "\n",
    "The diffusion equation can be used to represent a great deal of natural\n",
    "and environmental processes. It was introduced by Fourier in 1822 to\n",
    "calculate the distribution of the temperature in materials and has later\n",
    "been applied by Fick to material science. The mathematical expression\n",
    "that we will derive can be used to model e.g. heat transfer in the\n",
    "earth's crust, soil evolution, transport of contaminant in an aquifer or\n",
    "in the atmosphere, erosion of mountain ranges, the evolution of glaciers\n",
    "and many other phenomena. But before describing the equation directly,\n",
    "we will investigate what diffusion actually means.\n",
    "\n",
    "Note: Lecture notes on diffusion are partly based on Prof. Dr. Frédéric Herman's course on geophysical processes. \n",
    "@author: Benjamin Campforts\n",
    "\n",
    "#### Diffusion, what does it mean? \n",
    "The following movie illustrates Brownian motion, we can\n",
    "see that the equation to derive for diffusion must enable us to\n",
    "represent the movement of molecules from a high concentration zone to a\n",
    "zone of low concentration (Movie 1): \n",
    "\n",
    "[![Brownian motion](../../media/diffusion.png)](https://www.youtube.com/watch?v=UhL9OsRSKO8 \"Brownian motion\")\n",
    "\n",
    "**Movie 1:** Brownian motion causes food dye molecules to move throughout the water\n",
    "\n",
    "If we simplify this is a graph, we would get the following: \n",
    "<img src=\"../../media/Diff_Fig1.png\" style=\"width:3in;height:2in\" />\n",
    "\n",
    "**Figure 1**: Diffusion is the movement of molecules from a high\n",
    "concentration zone to a zone of low concentration due to random\n",
    "processes. C represents the concentration; X is the horizontal distance\n",
    "and q is the net particle flow.\n",
    "\n",
    "Due to diffusion, the particles move from the black zone to the grey\n",
    "zone. This can be explained by the fact that each particle can move at\n",
    "any moment in any direction, over a given distance. In one dimension, a\n",
    "particle can move to the left or to the right with equal probability,\n",
    "and this as well in the gray region as the black region. However, at the\n",
    "transition from the black zone to the grey zone, the probability of\n",
    "seeing particles move from left to right is much larger than the\n",
    "opposite (because there are much more black particles). This causes a\n",
    "particle transfer that depends on the difference of concentration $\\Delta C$ and\n",
    "the distance that the particle must travel $\\Delta x$, where $\\Delta C$ is the\n",
    "difference of concentration in a transition zone of length $\\Delta x$.\n",
    "Therefore, we can see that the flow of particles (i.e. the number of\n",
    "particles passing through per unit surface and time (in 2D, mol m<sup>-1</sup>\n",
    "s<sup>-1</sup>) will depend on the concentration gradient. Over time,\n",
    "the concentration changes as illustrated in Figure 2.\n",
    "\n",
    "<img src=\"../../media/Diff_Fig2.png\" style=\"width:3in;height:2in\" />\n",
    "\n",
    "**Figure 2:** Concentration changes over time due to diffusion\n",
    "\n",
    "### Exercise 1: modeling the random movement of particles\n",
    "\n",
    "Start by importing some libraries we will use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of the first exercise is to replicate the process of diffusion by modeling the random movement of particles.\n",
    "\n",
    "First, create a vector (called `xp`) that contains 100000 particles, having a random value between -20 and 20. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Represent this graphically using a histogram (use the `plt.hist()` function) and using the number of bins `xbins` (see below) in which you calculate the frequency (that is, the number of particles that is in each bin). Complete the code block below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xbins = np.arange(-200,200,2) # The number of bins in which you calculate the frequency (i.e. the number of particles that is in each bin.) \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, create a loop in which:\n",
    "- each particle moves over a random distance (positive or negative). This distance varies from -2 to 2. Use the `random` function. \n",
    "- Run the loop 100 times (`nbT` =100). You will need two for loops (also referred to as nested for loops). \n",
    "- Plot the resulting distribution every 200 iterations (using the code you coded above). Fix the limits of the y axes (`plt.ylim()`) to better understand how the distribution changes with time. \n",
    "- Extra: can you find a way to check how long it takes to execute the loop (a timer)? \n",
    "\n",
    "The code should have the following structure: \n",
    "~~~\n",
    "nbT = ...\n",
    "for t ...:\n",
    "    for i ...:\n",
    "        xp[i] ...\n",
    "    if t%50==0:\n",
    "        plt.figure()\n",
    "        plt.hist(xp, bins=xbins)\n",
    "        plt.title('Time is: ' +str(t))\n",
    "        plt.xlabel('Horizontal distance')\n",
    "        plt.ylabel('Number of particles')\n",
    "        plt.ylim((0, 10e3))\n",
    "        plt.show()\n",
    "\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you notice, this code is very slow. Vectorize your problem using numpy arrays to speed up the calculation and get rid of the inner for loop (you can keep the loop taking care of time). Change `nbT` to 10000 iterations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that you have an efficient solution, answer the following questions:\n",
    "\n",
    "1. Describe the evolution of the particles. How does the shape of the histogram evolve?\n",
    "2. Why does the evolution of the histogram over time slows down?\n",
    "3. What happens if the number of particles is reduced to 1000?\n",
    "4. What happens if the distance of the displacement now varies randomly between -10 and 10? Why is this the case?\n",
    "\n",
    "\n",
    "## Derivation of the diffusion equation\n",
    "\n",
    "In the previous exercise, we modeled a particle transfer assuming a random particle shift. The solution showed that the change in particle distribution depends on the concentration difference $\\Delta C$ and the distance the particle must travel $\\Delta x$, where $\\Delta C$ is the difference in concentration in the transition zone of length $\\Delta x$ (Figure 2).\n",
    "\n",
    "From these observations, we can thus conclude that particle flow, i.e. the number of particles passing through the side of an infinitesimal block per unit of time $\\mathrm{(mol \\: m^{-1} s^{-1})}$, will depend on the concentration gradient (Figure 2).\n",
    "\n",
    "We can therefore say that the flux, $q$ is defined by:\n",
    "\n",
    "$$q = -D\\frac{\\Delta C}{\\Delta x} \\label{eq:1} \\tag{1}$$\n",
    "\n",
    "\n",
    "where $D$ corresponds to the diffusion coefficient $\\mathrm{(m^{2} s^{-1})}$, or the diffusivity. *C* represents the concentration or the number of elements in a 2-dimensional infinitesimal block $\\mathrm{(mol \\: m^{-2})}$. The diffusion coefficient will vary from one problem to another and defines the speed of particle transfer. Now, we would also like to know how the *concentration* changes during the calculations. Let's take an infinitesimal block with an incoming flow, and an outgoing flow (Figure 4).\n",
    "\n",
    "<img src=\"../../media/Diff_Fig4.png\" style=\"width:3in;height:2in\"  />\n",
    "\n",
    "**Figure 4:** Infinitesimal block with an incoming flow, and an\n",
    "outgoing flow.\n",
    "\n",
    "As the concentration varies in the *X* direction, the flow will be\n",
    "different at the input and at the output of the block. The difference of\n",
    "the number of particles in the block can therefore be derived from the\n",
    "flux:\n",
    "\n",
    "$$ \\Delta (number \\: of \\: particles) = (q_x -q_{x+dx}\\Delta Y \\Delta t)\\label{eq:A1} \\tag{A1}$$\n",
    "\n",
    "Note that the term $\\Delta t$ appears because the flux dimensions are in $\\mathrm{(mol \\: m^{-1} s^{-1})}$ and that $\\Delta (number \\: of \\: particles)$ is in mol.\n",
    "\n",
    "We also know that the concentration change over an infinitesimally small\n",
    "unit of time corresponds to the change in the number of particles for a\n",
    "given volume:\n",
    "\n",
    "$$ \\Delta C=  \\frac{\\Delta (number \\: of \\: particles) }{\\Delta X \\Delta Y} \\label{eq:A2} \\tag{A2}$$\n",
    "\n",
    "which gives:\n",
    "\n",
    "$$ \\Delta (number \\: of \\: particles) = \\Delta C \\Delta X \\Delta Y \\label{eq:A3} \\tag{A3}$$\n",
    "\n",
    "By combining Equation (\\ref{eq:A1}) and (\\ref{eq:A3}), we can write:\n",
    "\n",
    "$$ (q_x -q_{x+dx}) \\Delta Y \\Delta t = \\Delta C \\Delta X \\Delta Y  \\label{eq:A4} \\tag{A4}$$\n",
    "\n",
    "$$ (q_x -q_{x+dx}) \\Delta t = \\Delta C \\Delta X  \\label{eq:A5} \\tag{A5}$$\n",
    "\n",
    "\n",
    "$$ \\frac{\\Delta C}{\\Delta t} = \\frac{q_x -q_{x+dx}}{\\Delta X} \\label{eq:A6} \\tag{A6}$$\n",
    "\n",
    "Using the definition of a differential equation:\n",
    "\n",
    "$$ \\frac{\\delta q}{\\delta x} = \\frac{(q_{x+dx} -q_x)}{\\Delta X} \\label{eq:A7} \\tag{A7}$$\n",
    "\n",
    "We obtain the following equation (note the use of the $\\partial$ symbol: we solve a PDE):\n",
    "\n",
    "$$ \\frac{\\partial C}{\\partial t} = -\\frac{\\partial q}{\\partial x}  \\label{eq:2} \\tag{2}$$\n",
    "By combining Eqs.(\\ref{eq:1}) and (\\ref{eq:2}), we finally obtain the heat\n",
    "equation: \n",
    "\n",
    "$$ \\frac{\\partial C}{\\partial t} = D\\frac{\\partial^2 C}{\\partial x^2}  \\label{eq:3} \\tag{3}$$\n",
    "\n",
    "which depends only on the curvature (i.e. the second derivative) of the concentration and the diffusion constant. Therefore, it is sufficient to know the diffusion coefficient $D$ (which can be measured) and to measure the curvature to estimate the change in concentration over time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2: Change in concentration due to diffusion\n",
    "\n",
    "You will now solve the diffusion equation and write a code that allows us to solve this equation. The change in concentration will be calculated over a distance $Lx$. There are two ways to do this, we can either calculate the second derivative directly (i.e. the curvature), or do it in two steps by calculating the flux (i.e. the first derivative of the concentration) and then the derivative of the flux. We will use the second method because it is easier to calculate a first derivative than a second derivative.\n",
    "\n",
    "To be able to do the calculation we also need an initial condition (i.e. the starting concentration) and boundary conditions (i.e. the concentration in $x = 0$ and $x = Lx$). Finally, you will have to choose a time step that is small enough. The [Von Neumann stability analysis](https://en.wikipedia.org/wiki/Von_Neumann_stability_analysis) prescribes that $\\Delta t$ must be smaller than $\\frac{\\Delta X^2}{2D}$. \n",
    "\n",
    "Now, try to solve the diffusion equation through discretization of ([Eq. 3](#mjx-eqn-eq:3)). \n",
    "Make the following assumptions: \n",
    "1. The initial condition:\n",
    "    -   $C(x<=\\frac{Lx}{2}) = 500 \\: \\mathrm{(mol \\: m^{-2})}$\n",
    "    -   $C(x>\\frac{Lx}{2}) = 0 \\: \\mathrm{(mol \\: m^{-2})}$\n",
    "    -   $Lx = 30  \\: \\mathrm{m}$ or $Lx = 300   \\: \\mathrm{m}$\n",
    "    -   $D  = 100 \\:  \\mathrm{(m^{2} s^{-1})}$\n",
    "    -   $\\Delta x = 0.1  \\: \\mathrm{m}$\n",
    "\n",
    "2. Assumtions regarding the boundary conditions:\n",
    "    -   $C(x=0)  = 500 \\: \\mathrm{(mol \\: m^{-2})}$\n",
    "    -   $C(x=Lx) = 0 \\: \\mathrm{(mol \\: m^{-2})}$\n",
    "    \n",
    "The code to solve this exercise must have the following structure:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "~~~\n",
    "#physics\n",
    "D = \n",
    "Lx = \n",
    "time =\n",
    "\n",
    "#numerical properties\n",
    "dx = \n",
    "x = np.arange(...)\n",
    "nx =\n",
    "nt = \n",
    "nout = \n",
    "\n",
    "# initial condition\n",
    "# Choose an initial condition where C = C1 when x <= (Lx/2) and C = 0 when x > (Lx/2)\n",
    "C1 = \n",
    "C2 = \n",
    "C  =\n",
    "C[x<=Lx/2] = \n",
    "C[x>Lx/2] = \n",
    "\n",
    "# Plot the initial concentration \n",
    "# plt.figure()\n",
    "# plt.plot(x,C)\n",
    "# plt.title('Initial condition')\n",
    "# plt.show()\n",
    "\n",
    "# impose a condition on the time step (Von Neumann stability criterion)\n",
    "dt = dx*dx/D/2.1\n",
    "print(dt)\n",
    "\n",
    "#model run: solve the heat equation and plot the result.\n",
    "\n",
    "# - make a time loop\n",
    "for t in range(0,nt):\n",
    "    # - in this loop, first calculate the flux by discretizing equation (1), \n",
    "    #   try to use vectorized code (eg. using numPy diff statement)\n",
    "    q = \n",
    "    \n",
    "    # - Update the new concentration (Eq. 2, without changing the boundary values)\n",
    "    #  Careful: which nodes do you have to update now? \n",
    "    C[...] = \n",
    "    \n",
    "    # - plot intermediate results, but only for every 100 time steps\n",
    "    if t%100==0:\n",
    "        #plt.figure()\n",
    "        plt.plot(x,C)\n",
    "        plt.title('Time is: ' +str(t))\n",
    "        \n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, try to answer the following questions: \n",
    "\n",
    "1. What is the shape of the concentration in equilibrium?\n",
    "2. How long does it take to reach equilibrium?\n",
    "3. Let *Lx* vary between 30 and 300, and *D* between 20 and 500. How does the time change to arrive at equilibrium according to *L* and *D*?\n",
    "\n",
    "To be able to answer the questions, you can modify your code to assume a condition on the concentration that defines when the solution will have reached a state of equilibrium. To implement this, replace the loop 'for' with a while loop: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "~~~\n",
    "Cp=C\n",
    "it=0\n",
    "diff=1e6\n",
    "\n",
    "while ...: #(diff > 1e-4)\n",
    "    it +=1 \n",
    "    #update the time\n",
    "    #calculate the flow with the discretized equation (eq. 1)\n",
    "    #calculate the new concentration (eq. 2)(without changing the BC's)\n",
    "\n",
    "    #check if the solution changes\n",
    "    diff = # sum of absolute difference between Cp and C\n",
    "    Cp = C\n",
    "    #plotting (only every 10000 iterations)\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Practice your skills: Eyjafjallajokull- Part 1\n",
    "\n",
    "<img src=\"../../media/Eyjafjallajokull.jpg\" />\n",
    "\n",
    "**Figure** \"On March 2010 Eyjafjallajokull volcano in Iceland exploded into life, spewing lava, magma, rock and clouds of ash into the sky above it. The disaster grounded airlines, stranding holidaymakers and business passengers across Europe and North America. While many could only watch the crisis unfolding in hotels and airports, photographer Ragnar Th. Sigurdsson chose to fly into the epicentre on a mission to record one of nature's most deadly phenomena\" The Telegraph (c); Picture: RAGNAR TH. SIGURDSSON / BARCROFT \n",
    "\n",
    "Exercise: 1-D  diffusion\n",
    "The Eyjafjallajökull volcano is located at the Southern Iceland coast and 2000 km from Brussels. Consider a one-dimensional case with a domain length of 5000 km. The volcano itself is situated at 2220 km from the left boundary of the simulation domain. Brussels is situated at 4220 km from the left boundary of the simulation domain. Choose a spatial resolution of 20 km. In the next couple of steps, you will calculate the time required to obtain a specific ash concentration above Brussels. \n",
    "\n",
    "<img src=\"../../media/Situation1D.png\" />\n",
    "\n",
    "**Figure:** 1D situation sketch\n",
    "\n",
    "Solve the spread of ash using the diffusion equation (Eq. 3)\n",
    "\n",
    "Define the model parameters: \n",
    "- set the diffusivity to 25 km$^2$/h\n",
    "- define the model domain: total length is 5000 km, spatial resolution is 20 km\n",
    "- calcualte the location (index in the np array) of the volcano and of Brussels and call the variables respectively `ind_vol` and `ind_Bru`\n",
    "- Set the initial conditions (C): at the volcano the concentration is 100 ppm, over the rest of the domain the concentration is 0 ppm.\n",
    "- The Eyjafjallajökull volcano produced ashes almost continuously during a couple of weeks. Start from the initial condition above but now add 100 ppm ashes per hour to the volcano grid cell as a source term. \n",
    "- Assume Dirichlet boundary conditions (0 ppm at 0 km and 0 ppm at 3000 km)\n",
    "- Plot the initial concentration, also indicate the location of Brussels on the plot (HINT use `plt.scatter()`)\n",
    "- Calculate and print out the time step (dt) using the CFL criterion\n",
    "\n",
    "The code must have the following structure: \n",
    "\n",
    "~~~\n",
    "#physics\n",
    "D = \n",
    "Lx = \n",
    "time =\n",
    "\n",
    "#numerical properties\n",
    "dx = \n",
    "x = \n",
    "nx =\n",
    "nt = \n",
    "nout = \n",
    "\n",
    "# Location of volcano and Brussels\n",
    "ind_vol= \n",
    "ind_Bru= \n",
    "\n",
    "C_ini  = \n",
    "C_rate = \n",
    "Cstart = 0\n",
    "Cend  = 0\n",
    "C =np.zeros(x.shape)\n",
    "\n",
    "C[0] = Cstart\n",
    "C[ind_vol] = C_ini\n",
    "C[-1] = Cend\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(x,C)\n",
    "plt.scatter(x[ind_Bru],C[ind_Bru],c='r')\n",
    "\n",
    "dt = dx*dx/D/2.5\n",
    "print('dt is: ' + str(dt) + 'hours')\n",
    "\n",
    "it = 0\n",
    "plt.figure()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- After how many hours do we get 5 ppm ash aerosols in Brussels?\n",
    "- Is this realistic? \n",
    "\n",
    "Use the solution derived before to solve this question using a while loop. Plot the output every 100 iterations. \n",
    "\n",
    "The code must have the following structure: \n",
    "\n",
    "~~~\n",
    "\n",
    "it =0 \n",
    "\n",
    "while ...:\n",
    "    it...\n",
    "    q =...\n",
    "    C...\n",
    "    \n",
    "    # Source term\n",
    "    C[ind_vol] += C_rate*dt\n",
    "    \n",
    "    # Boundary conditions\n",
    "     C[0] = Cstart\n",
    "     C[-1] = Cend\n",
    "\n",
    "    \n",
    "    if ...:        \n",
    "        plt.plot(x,C)\n",
    "        plt.scatter(x[ind_Bru],C[ind_Bru],c='r')\n",
    "        plt.title('Time is: ' + str(it*dt) + ' sec')    \n",
    "        plt.show()\n",
    "            \n",
    "print('Concentration reached after: ' + str(int(it*dt)) + ' hours')\n",
    "print('or : ' + str(int(it*dt/24)) + ' days')\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
